module LTK_universe
using LinearAlgebra
using SparseArrays
using Random
#export Universe


import ..Gaugefields:
    GaugeFields,
    GaugeFields_1d,
    IdentityGauges,
    RandomGauges,
    set_wing!,
    substitute!,
    calc_GaugeAction,
    SU3GaugeFields,
    SU3GaugeFields_1d,
    SU2GaugeFields,
    SU2GaugeFields_1d,
    SUNGaugeFields,
    SUNGaugeFields_1d,
    Oneinstanton,
    evaluate_wilson_loops!,
    U1GaugeFields,
    U1GaugeFields_1d,
    apply_smearing,
    calc_smearingU
import ..Gaugefields

import ..Fermionfields:
    FermionFields,
    WilsonFermion,
    StaggeredFermion,
    substitute_fermion!,
    gauss_distribution_fermi!,
    set_wing_fermi!
import ..Fermionfields
import ..Actions:
    GaugeActionParam,
    FermiActionParam,
    Setup_Gauge_action,
    Setup_Fermi_action,
    GaugeActionParam_standard,
    FermiActionParam_Wilson,
    show_parameters_action,
    FermiActionParam_WilsonClover,
    FermiActionParam_Staggered,
    GaugeActionParam_autogenerator,
    SmearingParam_multi,
    SmearingParam_single,
    Nosmearing,
    FermiActionParam_Domainwall
import ..LieAlgebrafields: LieAlgebraFields, clear!, add_gaugeforce!
import ..LieAlgebrafields
import ..Rand: Random_LCGs
import ..System_parameters: Params
import ..Diracoperators:
    DdagD_operator, Domainwall_operator, make_densematrix, D5DW_Domainwall_operator
import ..Diracoperators
import ..Wilsonloops:
    make_loopforactions, Wilson_loop_set, make_originalactions_fromloops, make_cloverloops
import ..Verbose_print: Verbose_level, Verbose_3, Verbose_2, Verbose_1
import ..IOmodule: loadU
import ..SUN_generator: Generator

import ..RationalApprox: calc_exactvalue, calc_Anϕ
import ..ILDG_format: ILDG, load_gaugefield, load_gaugefield!, save_binarydata
import ..Bridge_format: load_BridgeText!
import ..Othermethods: tdlogdet
import ..CGmethods: bicg, cg, shiftedcg
import ..Rhmc: get_order, get_β, get_α, get_α0, get_β_inverse, get_α_inverse, get_α0_inverse
#import ..Smearing:apply_smearing



"""
Your universe is described in this type.
"""
struct Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp}
    NX::Int64
    NY::Int64
    NZ::Int64
    NT::Int64
    NV::Int64
    NC::Int64
    Nwing::Int64
    Dirac_operator::Union{Nothing,String}
    U::Array{Gauge,1}
    Uold::Array{Gauge,1}
    p::Array{Lie,1}
    φ::Union{Nothing,Fermi}
    ξ::Union{Nothing,Fermi}
    η::Union{Nothing,Fermi}
    gparam::GaugeP
    fparam::Union{Nothing,FermiP}
    BoundaryCondition::Array{Int8,1}
    initial::String
    quench::Bool
    NDFALG::Int64





    _temporal_gauge::Array{Gauge_temp,1}
    _temporal_fermi::Union{Nothing,Array{Fermi,1}}
    _temporal_algebra::Array{Lie,1}
    ranf::Random_LCGs
    verboselevel::Int8
    kind_of_verboselevel::Verbose_level
    #rng::MersenneTwister

end





function set_β!(univ::Universe, β)
    if typeof(univ.gparam) == GaugeActionParam_standard
        univ.gparam.β = β
        #univ.gparam = GaugeActionParam_standard(β,univ.gparam.NTRACE)
    elseif typeof(univ.gparam) == GaugeActionParam_autogenerator
        #univ.gparam =GaugeActionParam_autogenerator(univ.gparam.βs * (univ.gparam.β/β),β,univ.gparam.numactions,univ.gparam.NTRACE,univ.gparam.loops,univ.gparam.staples)
        univ.gparam.β = β
        univ.gparam.βs .= 0#univ.gparam.βs * (univ.gparam.β/β)
        univ.gparam.βs[1] = β
    else
        error("$(typeof(univ.gparam)) is not supported!")
    end
end

function set_β!(univ::Universe, βs::Array{T,1}) where {T<:Number}
    if typeof(univ.gparam) == GaugeActionParam_standard
        univ.gparam.β = βs[1]
        #univ.gparam = GaugeActionParam_standard(β,univ.gparam.NTRACE)
    elseif typeof(univ.gparam) == GaugeActionParam_autogenerator
        @assert length(univ.gparam.βs) == length(βs) "univ.gparam.βs = $(univ.gparam.βs) and βs = $βs"
        #univ.gparam =GaugeActionParam_autogenerator(univ.gparam.βs * (univ.gparam.β/β),β,univ.gparam.numactions,univ.gparam.NTRACE,univ.gparam.loops,univ.gparam.staples)
        univ.gparam.β = βs[1]
        univ.gparam.βs .= βs
        #univ.gparam.βs = univ.gparam.βs * (univ.gparam.β/β)
    else
        error("$(typeof(univ.gparam)) is not supported!")
    end
end

function get_β(univ::Universe)
    if typeof(univ.gparam) == GaugeActionParam_standard
        return univ.gparam.β
    elseif typeof(univ.gparam) == GaugeActionParam_autogenerator
        return univ.gparam.β
    else
        error("$(typeof(univ.gparam)) is not supported!")
    end
end

function set_βs!(univ::Universe, βs)
    if typeof(univ.gparam) == GaugeActionParam_standard
        univ.gparam.β = βs[1]
    elseif typeof(univ.gparam) == GaugeActionParam_autogenerator
        univ.gparam.β = βs[1]
        univ.gparam.βs = βs[:]
    else
        error("$(typeof(univ.gparam)) is not supported!")
    end
end


include("default.jl")

function Universe()
    file = "default.jl"
    return Universe(file)
end

function Universe(p::Params)
    L = p.L
    β = p.β
    NTRACE = p.NC
    #rng = MersenneTwister(p.randomseed, (0, 1002, 0, 1))

    if p.use_autogeneratedstaples
        #loops = make_loopforactions(p.couplinglist,L)

        #println(loops)
        loops = make_originalactions_fromloops(p.coupling_loops)

        gparam =
            GaugeActionParam_autogenerator(p.couplingcoeff, loops, p.NC, p.couplinglist)
    else
        gparam = GaugeActionParam_standard(β, NTRACE)
    end

    if p.Dirac_operator == nothing
        fparam = nothing
    else
        if p.Dirac_operator == "Wilson"
            if p.smearing_for_fermion == "nothing"
                fparam = FermiActionParam_Wilson(
                    p.hop,
                    p.r,
                    p.eps,
                    p.Dirac_operator,
                    p.MaxCGstep,
                    p.quench,
                )
            else
                fparam = FermiActionParam_Wilson(
                    p.hop,
                    p.r,
                    p.eps,
                    p.Dirac_operator,
                    p.MaxCGstep,
                    p.quench,
                    smearingparameters = "stout",
                    loops_list = p.stout_loops,
                    coefficients = p.stout_ρ,
                    numlayers = p.stout_numlayers,
                    L = p.L,
                )
            end
        elseif p.Dirac_operator == "WilsonClover"
            #if p.NC == 2
            #    error("You use NC = 2. But WilsonClover Fermion is not supported in SU2 gauge theory yet. Use NC = 3.")
            #end
            NV = prod(p.L)
            #CloverFμν = zeros(ComplexF64,p.NC,p.NC,NV,6)
            inn_table = zeros(Int64, NV, 4, 2)
            internal_flags = zeros(Bool, 2)
            _ftmp_vectors = Array{Array{ComplexF64,3},1}(undef, 6)
            for i = 1:6
                _ftmp_vectors[i] = zeros(ComplexF64, p.NC, NV, 4)
            end

            _is1 = zeros(Int64, NV)
            _is2 = zeros(Int64, NV)

            if p.NC ≥ 4
                SUNgenerator = Generator(p.NC)
            else
                SUNgenerator = nothing
            end

            _cloverloops = Array{Wilson_loop_set,2}(undef, 3, 4)
            for μ = 1:3
                for ν = μ+1:4
                    _cloverloops[μ, ν] = make_cloverloops(μ, ν)
                end
            end


            #fparam = FermiActionParam_WilsonClover(p.hop,p.r,p.eps,p.Dirac_operator,p.MaxCGstep,p.Clover_coefficient,CloverFμν,
            #                internal_flags,inn_table,_ftmp_vectors,_is1,_is2,
            #                p.quench)
            if p.smearing_for_fermion == "nothing"
                fparam = FermiActionParam_WilsonClover(
                    p.hop,
                    p.r,
                    p.eps,
                    p.Dirac_operator,
                    p.MaxCGstep,
                    p.Clover_coefficient,
                    internal_flags,
                    inn_table,
                    _ftmp_vectors,
                    _is1,
                    _is2,
                    p.quench,
                    SUNgenerator,
                    _cloverloops,
                )
            else
                error("stout for Wilson with clover is not supported yet!")
                fparam = FermiActionParam_WilsonClover(
                    p.hop,
                    p.r,
                    p.eps,
                    p.Dirac_operator,
                    p.MaxCGstep,
                    p.Clover_coefficient,
                    internal_flags,
                    inn_table,
                    _ftmp_vectors,
                    _is1,
                    _is2,
                    p.quench,
                    SUNgenerator,
                    _cloverloops,
                    smearingparameters = "stout",
                    loops_list = p.stout_loops,
                    coefficients = p.stout_ρ,
                    numlayers = p.stout_numlayers,
                    L = p.L,
                )
            end
        elseif p.Dirac_operator == "Staggered"
            if p.smearing_for_fermion == "nothing"
                fparam = FermiActionParam_Staggered(
                    p.mass,
                    p.eps,
                    p.Dirac_operator,
                    p.MaxCGstep,
                    p.quench,
                    p.Nf,
                )
            else
                fparam = FermiActionParam_Staggered(
                    p.mass,
                    p.eps,
                    p.Dirac_operator,
                    p.MaxCGstep,
                    p.quench,
                    p.Nf,
                    smearingparameters = "stout",
                    loops_list = p.stout_loops,
                    coefficients = p.stout_ρ,
                    numlayers = p.stout_numlayers,
                    L = p.L,
                )
            end
        elseif p.Dirac_operator == "Domainwall"
            if p.smearing_for_fermion == "nothing"
                fparam = FermiActionParam_Domainwall(
                    p.Domainwall_N5,
                    p.r,
                    p.Domainwall_M,
                    p.Domainwall_m,
                    p.Domainwall_ωs,
                    p.Domainwall_b,
                    p.Domainwall_c,
                    p.eps,
                    p.Dirac_operator,
                    p.MaxCGstep,
                    p.quench,
                )
            else
                fparam = FermiActionParam_Domainwall(
                    p.Domainwall_N5,
                    p.r,
                    p.Domainwall_M,
                    p.Domainwall_m,
                    p.Domainwall_ωs,
                    p.Domainwall_b,
                    p.Domainwall_c,
                    p.eps,
                    p.Dirac_operator,
                    p.MaxCGstep,
                    p.quench,
                    smearingparameters = "stout",
                    loops_list = p.stout_loops,
                    coefficients = p.stout_ρ,
                    numlayers = p.stout_numlayers,
                    L = p.L,
                )
            end
            #error("Not implemented! here is in LTK_universe.jl")
        else
            error(p.Dirac_operator, " is not supported!")
        end
    end

    univ = Universe(
        L,
        gparam,
        p.Nwing,
        fparam,
        p.BoundaryCondition,
        p.initial,
        p.NC,
        p.verboselevel,
        p.load_fp,
        p.loadU_format,
        p.julian_random_number,
    )

end




"""
```Universe(file)```
- file: file name of the input file.

Make your universe. The input file is loaded.

Undefined parameters in your input file are defined with the default values.

The default values are as follows.

```julia
    L = (4,4,4,4)
    β = 6
    NTRACE = 3
    #gparam = Setup_Gauge_action(β)
    gparam =  GaugeActionParam_standard(β,NTRACE)

    BoundaryCondition=[1,1,1,-1]
    Nwing = 1
    initial="cold"
    NC =3


    hop= 0.141139 #Hopping parameter
    r= 1 #Wilson term
    eps= 1e-19
    Dirac_operator= "Wilson"
    MaxCGstep= 3000

    fparam = FermiActionParam_Wilson(hop,r,eps,Dirac_operator,MaxCGstep)

```
"""
function Universe(file)
    if file != "default.jl"
        include(pwd() * "/" * file)
    end
    univ = Universe(L, gparam, Nwing, fparam, BoundaryCondition, initial, NC)

end

function show_parameters(univ::Universe)
    println("""
        L = $((univ.NX,univ.NY,univ.NZ,univ.NT))
        """)

    show_parameters_action(univ.gparam)

    println("""
        BoundaryCondition= $(univ.BoundaryCondition)
        Nwing = $(univ.Nwing)
        initial= $(univ.initial)
        NC =$(univ.NC)
        """)

    if univ.fparam != nothing
        show_parameters_action(univ.fparam)
    end
end

"""
```Universe(L::Tuple,gparam::GaugeActionParam,initial="cold",fparam=nothing)```

- L: system size (NX,NY,NZ,NT)
- gparam: parameters for gauge actions
- [initial]: initial Gauge configuration
- [fparam]: parameters for fermion actions
"""
function Universe(L::Tuple, gparam::GaugeActionParam, initial = "cold", fparam = nothing)
    Nwing = 1
    BoundaryCondition = [1, 1, 1, -1]
    NC = 3
    Universe(L, gparam, Nwing, fparam, BoundaryCondition, initial, NC)
end


function Universe(
    L,
    gparam,
    Nwing,
    fparam,
    BoundaryCondition,
    initial,
    NC,
    verboselevel,
    load_fp,
    loadU_format,
    julian_random_number,
)
    #(L::Tuple,gparam::GaugeActionParam;
    #   Nwing = 1,fparam=nothing,
    #  BoundaryCondition=[1,1,1,-1],initial="cold",NC =3)

    if verboselevel == 1
        kind_of_verboselevel = Verbose_1(load_fp)
    elseif verboselevel == 2
        kind_of_verboselevel = Verbose_2(load_fp)
    elseif verboselevel == 3
        kind_of_verboselevel = Verbose_3(load_fp)
    end


    ranf = Random_LCGs(1, julian_random_number = julian_random_number)

    NX = L[1]
    NY = L[2]
    NZ = L[3]
    NT = L[4]
    NV = NX * NY * NZ * NT
    NDFALG = ifelse(NC == 1, 1, NC^2 - 1)

    num_tempfield_g = 5
    num_tempfield_f = 4

    if fparam != nothing && fparam.Dirac_operator == "WilsonClover"
        numbasis = NC^2 - 1
        num_tempfield_g += 4 + 4 + 2numbasis
    elseif fparam != nothing # && fparam.Dirac_operator == "Staggered"
        num_tempfield_f += 4
    end

    if fparam != nothing && fparam.Dirac_operator == "Staggered"
        if fparam.Nf != 4 && fparam.Nf != 8
            N_action = get_order(fparam.rhmc_action)
            N_MD = get_order(fparam.rhmc_MD)
            num_tempfield_f += maximum((N_action, N_MD)) + 1
        end
    end

    if fparam != nothing && fparam.smearing != nothing
        num_tempfield_g += 4
    end

    if typeof(gparam) == GaugeActionParam_autogenerator
        num_tempfield_g += 1
    end


    if NC == 3
        U = Array{SU3GaugeFields,1}(undef, 4)
        _temporal_gauge = Array{SU3GaugeFields_1d,1}(undef, num_tempfield_g)
    elseif NC == 2
        U = Array{SU2GaugeFields,1}(undef, 4)
        _temporal_gauge = Array{SU2GaugeFields_1d,1}(undef, num_tempfield_g)
    elseif NC ≥ 4
        U = Array{SUNGaugeFields{NC},1}(undef, 4)
        _temporal_gauge = Array{SUNGaugeFields_1d{NC},1}(undef, num_tempfield_g)
    elseif NC == 1
        U = Array{U1GaugeFields,1}(undef, 4)
        _temporal_gauge = Array{U1GaugeFields_1d,1}(undef, num_tempfield_g)
    end

    #        Uold = Array{GaugeFields,1}(undef,4)

    p = Array{LieAlgebraFields,1}(undef, 4)
    for μ = 1:4
        p[μ] = LieAlgebraFields(NC, NX, NY, NZ, NT)
    end

    if initial == "cold"
        println(".....  Cold start")
        for μ = 1:4
            U[μ] = IdentityGauges(NC, NX, NY, NZ, NT, Nwing)
        end
    elseif initial == "hot"
        println(".....  Hot start")
        for μ = 1:4
            U[μ] = RandomGauges(NC, NX, NY, NZ, NT, Nwing)
        end
    elseif initial == "one instanton"
        @assert NC == 2 "From one instanton start, NC should be 2!"
        U = Oneinstanton(NC, NX, NY, NZ, NT, Nwing)

    else #if initial == "file"
        println(".....  File start")
        println("File name is $initial")
        if loadU_format == "JLD"
            U = loadU(initial, NX, NY, NZ, NT, NC)
        elseif loadU_format == "ILDG"
            ildg = ILDG(initial)
            i = 1
            for μ = 1:4
                U[μ] = IdentityGauges(NC, NX, NY, NZ, NT, Nwing)
            end
            load_gaugefield!(U, i, ildg, L, NC)
        elseif loadU_format == "BridgeText"
            for μ = 1:4
                U[μ] = IdentityGauges(NC, NX, NY, NZ, NT, Nwing)
            end
            load_BridgeText!(initial, U, L, NC)
        else
            error("loadU_format should be JLD or ILDG")
        end

        #error("not supported yet.")
    end

    set_wing!(U)
    Uold = similar(U)
    substitute!(Uold, U)


    #=
    for μ=1:4
        set_wing!(U[μ])
        Uold[μ] = similar(U[μ])
        substitute!(Uold[μ],U[μ])
    end
    =#


    #_temporal_gauge = Array{GaugeFields,1}(undef,4)
    for i = 1:length(_temporal_gauge)
        _temporal_gauge[i] = GaugeFields_1d(NC, NX, NY, NZ, NT) #similar(U[1])
        #_temporal_gauge[i] = similar(U[1])
    end

    _temporal_algebra = Array{LieAlgebraFields,1}(undef, 1)
    for i = 1:length(_temporal_algebra)
        _temporal_algebra[i] = similar(p[1])
    end


    if fparam == nothing
        Dirac_operator = nothing
        φ = nothing
        η = nothing
        ξ = nothing
        _temporal_fermi = nothing
        quench = true
    else
        Dirac_operator = fparam.Dirac_operator
        φ = FermionFields(NC, NX, NY, NZ, NT, fparam, BoundaryCondition)
        η = similar(φ)
        ξ = similar(φ)
        _temporal_fermi = Array{FermionFields,1}(undef, num_tempfield_f)
        for i = 1:length(_temporal_fermi)
            _temporal_fermi[i] = similar(φ)
        end
        quench = fparam.quench

        if Dirac_operator == "Domainwall"
            gauss_distribution_fermi!(ξ)
            #A = D5DW_Domainwall_operator(U,ξ,fparam)
            A = Domainwall_operator(U, ξ, fparam)
            A_dense_PV = make_densematrix(A.D5DW_PV)
            A_dense = make_densematrix(A.D5DW)
            M = A_dense * inv(A_dense_PV)
            #e,v = eigen(A_dense)
            e, v = eigen(M)
            fp = open("eigenvalues_M.txt", "w")
            for ene in e
                println(fp, real(ene), "\t", imag(ene))
            end
            close(fp)
            #println(e)
            error("stop!")
            bicg(η, A, ξ)
            error("stop!")
        end

    end



    Gauge = eltype(U)
    Fermi = typeof(φ)
    Lie = eltype(p)
    GaugeP = typeof(gparam)
    FermiP = typeof(fparam)
    Gauge_temp = eltype(_temporal_gauge)



    return Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp}(
        NX,
        NY,
        NZ,
        NT,
        NV,
        NC,
        Nwing,
        Dirac_operator,
        U,
        Uold,
        p,
        φ,
        ξ,
        η,
        gparam,
        fparam,
        BoundaryCondition,
        initial,
        quench,
        NDFALG,
        _temporal_gauge,
        _temporal_fermi,
        _temporal_algebra,
        ranf,
        verboselevel,
        kind_of_verboselevel,
    )


end

function gauss_distribution(univ::Universe, nv)
    variance = 1
    nvh = div(nv, 2)
    granf = zeros(Float64, nv)
    for i = 1:nvh
        rho = sqrt(-2 * log(univ.ranf()) * variance)
        theta = 2pi * univ.ranf()
        granf[i] = rho * cos(theta)
        granf[i+nvh] = rho * sin(theta)
    end
    if 2 * nvh == nv
        return granf
    end

    granf[nv] = sqrt(-2 * log(univ.ranf()) * variance) * cos(2pi * univ.ranf())
    return granf
end






function expF_U!(
    U::Array{T,1},
    F::Array{N,1},
    Δτ,
    univ::Universe,
) where {T<:GaugeFields,N<:LieAlgebraFields}
    LieAlgebrafields.expF_U!(U, F, Δτ, univ._temporal_gauge, univ._temporal_algebra[1])
end

function calc_gaugeforce!(F::Array{N,1}, univ::Universe) where {N<:LieAlgebraFields}
    clear!(F)
    calc_gaugeforce!(F, univ.U, univ)
    #add_gaugeforce!(F,univ.U,univ._temporal_gauge,univ._temporal_algebra[1]) 
    return
end

function calc_gaugeforce!(
    F::Array{N,1},
    U::Array{T,1},
    univ::Universe,
) where {N<:LieAlgebraFields,T<:GaugeFields}
    clear!(F)
    calc_gaugeforce!(F, U, univ, univ.gparam)
    #add_gaugeforce!(F,U,univ._temporal_gauge,univ._temporal_algebra[1],univ.gparam) 
    return
end

function calc_gaugeforce!(
    F::Array{N,1},
    U::Array{T,1},
    univ::Universe,
    gparam::GP,
) where {N<:LieAlgebraFields,T<:GaugeFields,GP<:GaugeActionParam_autogenerator}
    clear!(F)
    add_gaugeforce!(F, U, univ._temporal_gauge, univ._temporal_algebra[1], gparam)
    return
end

function calc_gaugeforce!(
    F::Array{N,1},
    U::Array{T,1},
    univ::Universe,
    gparam::GP,
) where {N<:LieAlgebraFields,T<:GaugeFields,GP}
    clear!(F)
    add_gaugeforce!(F, U, univ._temporal_gauge, univ._temporal_algebra[1])
    return
end

function Gaugefields.calc_GaugeAction(univ::Universe)
    Sg, plaq = calc_GaugeAction(univ.U, univ.gparam, univ._temporal_gauge)
    return real(Sg), real(plaq)
end

function calc_Action(univ::Universe)
    Sg, plaq = calc_GaugeAction(univ.U, univ.gparam, univ._temporal_gauge)

    SP = univ.p * univ.p / 2
    #println("ek = $SP")
    S = Sg + SP
    #println("eym = $plaq")
    #println("e = $S")

    if univ.quench == false
        #if univ.fparam != nothing
        Sf = univ.η * univ.η
        S += Sf
    end
    return real(S), real(plaq)
end


function calc_IntegratedFermionAction(
    univ::Universe;
    debug = false,
    M = 16,
    m = 100,
    nc = 20,
    nonc = false,
)
    println("Making W^+W matrix...")
    NX = univ.U[1].NX
    NY = univ.U[1].NY
    NZ = univ.U[1].NZ
    NT = univ.U[1].NT
    NC = univ.U[1].NC
    T = eltype(univ._temporal_fermi)
    if T == StaggeredFermion
        NG = 1
    else
        NG = 4
    end
    Nsize = NX * NY * NZ * NT * NC * NG

    WdagW = zeros(ComplexF64, Nsize, Nsize)
    @time make_WdagWmatrix!(WdagW, univ)
    #@time WdagW = make_WdagWmatrix(univ)
    Sfnew = calc_IntegratedFermionAction(
        univ,
        WdagW,
        debug = debug,
        M = M,
        m = m,
        nc = nc,
        nonc = nonc,
    )
    return Sfnew
end


function calc_IntegratedFermionAction(
    univ::Universe,
    WdagW;
    debug = false,
    M = 16,
    m = 100,
    nc = 20,
    nonc = false,
)
    if debug

        #WdagW = DdagD_operator(univ.U,univ.η,univ.fparam)
        #println("Calculating approximated logdet")
        #@time Sfnew = -tdlogdet(WdagW,M,m,nc=nc,nonc=nonc)
        #exit()

        WdagWsp = sparse(WdagW)
        N, _ = size(WdagW)
        tempvecs = Array{Array{ComplexF64,1}}(undef, 3)
        for i = 1:3
            tempvecs[i] = zeros(ComplexF64, N)
        end
        println("Calculating approximated logdet")

        #@time Sfnew = -tdlogdet(WdagWsp,M2,m2,tempvecs,nc=nc2,nonc=nonc)

        #M2 = 100
        #m2 = 100
        #nc2 = 100
        @time Sfnew = -tdlogdet(WdagWsp, M, m, tempvecs, nc = nc, nonc = nonc)

        #exit()
        #=


        M = 8*2
        #M = 4#8*2
        m = 100
        NX = univ.U[1].NX
        #m = 20
        N,_ = size(WdagW)
        tempvecs = Array{Array{ComplexF64,1}}(undef,3)
        for i=1:3
            tempvecs[i] = zeros(ComplexF64,N)
        end

        nc = 20
        nonc = false

        @time ldet = logdet(WdagW)
        if nonc
            fp = open("deltavalue_m$(m)_nonc_NX$(NX)_6.dat","w")
        else
            fp = open("deltavalue_m$(m)_nc$(nc)_NX$(NX)_6.dat","w")
        end

        for M in [4,8,16,32,64,128]
            if nonc
                filename="delta_$(NX)$(NX)$(NX)$(NX)_$(M)_m$(m)_nonc_6.dat"
            else
                filename="delta_$(NX)$(NX)$(NX)$(NX)_$(M)_m$(m)_nc$(nc)_6.dat"
            end
            a = tdlogdet(WdagWsp,M,m,tempvecs,filename=filename,nc=nc,nonc=nonc)
            println("$M $a $ldet")
            println(fp,"$M $a $ldet")
        end
        close(fp)

        exit()
        @time Sfnew = -tdlogdet(WdagWsp,M,m,tempvecs)
        =#
    else
        println("Calculating logdet")
        @time Sfnew = -real(logdet(WdagW))
    end



    if univ.Dirac_operator == "Staggered"
        #if univ.fparam.Nf == 4
        #    Sfnew /= 2
        #end
        Sfnew /= (8 / univ.fparam.Nf)
    end
    return Sfnew
end

struct Wilsonloops_actions{Gauge_temp}
    loops::Array{Wilson_loop_set,1}
    loopaction::Gauge_temp
    temps::Array{Gauge_temp,1}
    numloops::Int64
    β::Float64
    NC::Float64
    couplinglist::Array{String,1}

    function Wilsonloops_actions(univ)
        couplinglist = ["plaq", "rect", "polyx", "polyy", "polyz", "polyt"]
        L = (univ.U[1].NX, univ.U[1].NY, univ.U[1].NZ, univ.U[1].NT)
        loops = make_loopforactions(couplinglist, L)
        loopaction = similar(univ._temporal_gauge[1])
        sutype = typeof(loopaction)
        temps = Array{sutype,1}(undef, 3)
        for i = 1:3
            temps[i] = similar(loopaction)
        end
        numloops = length(couplinglist)
        β = univ.gparam.β
        NC = univ.NC

        return new{sutype}(loops, loopaction, temps, numloops, β, NC, couplinglist)
    end
end


function calc_looptrvalues(w::Wilsonloops_actions, univ)
    numloops = w.numloops
    trs = zeros(ComplexF64, numloops)
    for i = 1:numloops
        evaluate_wilson_loops!(w.loopaction, w.loops[i], univ.U, w.temps[1:3])
        trs[i] = tr(w.loopaction)
    end
    return trs
end

function calc_looptrvalues_site(w::Wilsonloops_actions, univ)
    return
    numloops = w.numloops
    NX = univ.U[1].NX
    NY = univ.U[1].NY
    NZ = univ.U[1].NZ
    NT = univ.U[1].NT
    NC = univ.NC
    V = zeros(ComplexF64, NC, NC)

    trs = zeros(ComplexF64, numloops, NX, NY, NZ, NT)
    for it = 1:NT
        for iz = 1:NZ
            for iy = 1:NY
                for ix = 1:NX
                    for i = 1:numloops
                        evaluate_wilson_loops!(V, w.loops[i], univ.U, ix, iy, iz, it)
                        trs[i, ix, iy, iz, it] = tr(V)
                    end
                end
            end
        end
    end

    return trs
end

function get_looptrvalues(w::Wilsonloops_actions)
    return w.trs
end

function calc_IntegratedSf(w::Wilsonloops_actions, univ)
    Sf = calc_IntegratedFermionAction(univ)
    return Sf
end

function calc_gaugeSg(w::Wilsonloops_actions, univ)
    trs = calc_looptrvalues(w, univ)
    Sg = (-w.β / w.NC) * trs[1]
    return real(Sg), trs
end

function calc_trainingdata(w::Wilsonloops_actions, univ)
    Sf = calc_IntegratedSf(w, univ)
    Sg, trs = calc_gaugeSg(w, univ)
    return trs, Sg, Sf
end

function make_WdagWmatrix(univ::Universe)
    return make_WdagWmatrix(univ.U, univ._temporal_fermi, univ.fparam)
end

function make_WdagWmatrix!(WdagW, univ::Universe)
    make_WdagWmatrix!(WdagW, univ.U, univ._temporal_fermi, univ.fparam)
    return
end

function make_WdagWmatrix(univ::Universe, U::Array{T,1}) where {T<:GaugeFields}
    return make_WdagWmatrix(U, univ._temporal_fermi, univ.fparam)
end

function make_WdagWmatrix(
    U::Array{G,1},
    temps::Array{T,1},
    fparam,
) where {G<:GaugeFields,T<:FermionFields}
    x0 = temps[7]
    xi = temps[8]
    return make_WdagWmatrix(U, x0, xi, fparam)
end

function make_WdagWmatrix!(
    WdagW,
    U::Array{G,1},
    temps::Array{T,1},
    fparam,
) where {G<:GaugeFields,T<:FermionFields}
    x0 = temps[7]
    xi = temps[8]
    make_WdagWmatrix!(WdagW, U, x0, xi, fparam)
    return
end

function make_WdagWmatrix(
    U::Array{G,1},
    x0::T,
    xi::T,
    fparam,
) where {G<:GaugeFields,T<:FermionFields}
    NX = x0.NX
    NY = x0.NY
    NZ = x0.NZ
    NT = x0.NT
    NC = x0.NC
    if T == StaggeredFermion
        NG = 1
    else
        NG = 4
    end
    Nsize = NX * NY * NZ * NT * NC * NG

    WdagW = zeros(ComplexF64, Nsize, Nsize)
    make_WdagWmatrix!(WdagW, U, x0, xi, fparam)
    return WdagW
end

function make_WdagWmatrix!(
    WdagW,
    U::Array{G,1},
    x0::T,
    xi::T,
    fparam,
) where {G<:GaugeFields,T<:FermionFields}

    Fermionfields.clear!(x0)
    NX = x0.NX
    NY = x0.NY
    NZ = x0.NZ
    NT = x0.NT
    NC = x0.NC
    if T == StaggeredFermion
        NG = 1
    else
        NG = 4
    end
    Nsize = NX * NY * NZ * NT * NC * NG
    #WdagW = zeros(ComplexF64,Nsize,Nsize)
    WdagWoperator = DdagD_operator(U, x0, fparam)


    j = 0
    for α = 1:NG
        for it = 1:NT
            for iz = 1:NZ
                for iy = 1:NY
                    for ix = 1:NX
                        for ic = 1:NC
                            Fermionfields.clear!(x0)
                            j += 1
                            x0[ic, ix, iy, iz, it, α] = 1
                            #set_wing_fermi!(x0)

                            #WdagWx!(xi,U,x0,temps,fparam)
                            #println(xi*xi)

                            #Wx!(xi,U,x0,temps,fparam)
                            #println(xi*xi)

                            #Wx!(xi,U,x0,temps,fparam,(ix,iy,iz,it,α))
                            #println(xi*xi)

                            mul!(xi, WdagWoperator, x0, (ix, iy, iz, it, α))
                            #WdagWx!(xi,U,x0,temps,fparam,(ix,iy,iz,it,α))
                            #println("fast ",xi*xi)
                            #exit()
                            #display(xi)
                            #exit()
                            substitute_fermion!(WdagW, j, xi)

                        end
                    end
                end
            end
        end
    end


    #return WdagW

end


function construct_fermion_gauss_distribution!(
    univ::Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp},
) where {Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp}
    gauss_distribution_fermi!(univ.η, univ.ranf)
end

function construct_fermion_gauss_distribution!(
    univ::Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp},
) where {Gauge,Lie,Fermi<:StaggeredFermion,GaugeP,FermiP,Gauge_temp}
    gauss_distribution_fermi!(univ.η, univ.ranf)
    if univ.fparam.Nf == 4
        U, _... = calc_smearingU(univ.U, univ.fparam.smearing)

        evensite = false
        W = Diracoperators.Dirac_operator(U, univ.η, univ.fparam)
        mul!(univ.φ, W', univ.η)

        #Wdagx!(univ.φ,univ.U,univ.η,univ._temporal_fermi,univ.fparam)
        Fermionfields.clear!(univ.φ, evensite)

        bicg(univ.η, W', univ.φ, eps = univ.fparam.eps, maxsteps = univ.fparam.MaxCGstep)
    end
end

function construct_fermionfield_φ!(
    univ::Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp},
) where {Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp}
    U, _... = calc_smearingU(univ.U, univ.fparam.smearing)


    W = Diracoperators.Dirac_operator(U, univ.η, univ.fparam)
    mul!(univ.φ, W', univ.η)
    set_wing_fermi!(univ.φ)
end

function construct_fermionfield_φ!(
    univ::Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp},
) where {Gauge,Lie,Fermi<:StaggeredFermion,GaugeP,FermiP,Gauge_temp}
    U, _... = calc_smearingU(univ.U, univ.fparam.smearing)

    if univ.fparam.Nf == 4 || univ.fparam.Nf == 8

        W = Diracoperators.Dirac_operator(U, univ.η, univ.fparam)
        mul!(univ.φ, W', univ.η)
    else
        WdagW = Diracoperators.DdagD_operator(U, univ.η, univ.fparam)
        N = get_order(univ.fparam.rhmc_action)

        x = univ.φ
        vec_x = univ._temporal_fermi[end-N+1:end]
        for j = 1:N
            Fermionfields.clear!(vec_x[j])
        end

        # eta =  (MdagM)^{-alpha/2} phi -> phi = (MdagM)^{alpha/2} eta
        #(MdagM)^{alpha/2} eta ~ α0 eta + sum_k αk (MdagM + βk)^{-1} eta
        vec_β = get_β(univ.fparam.rhmc_action)
        vec_α = get_α(univ.fparam.rhmc_action)
        α0 = get_α0(univ.fparam.rhmc_action)
        shiftedcg(
            vec_x,
            vec_β,
            x,
            WdagW,
            univ.η,
            eps = univ.fparam.eps,
            maxsteps = univ.fparam.MaxCGstep,
        )
        Fermionfields.clear!(univ.φ)
        Fermionfields.add!(univ.φ, α0, univ.η)
        for j = 1:N
            αk = vec_α[j]
            Fermionfields.add!(univ.φ, αk, vec_x[j])
        end


    end

    set_wing_fermi!(univ.φ)
end


function construct_fermionfield_η!(
    univ::Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp},
) where {Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp}
    U, _... = calc_smearingU(univ.U, univ.fparam.smearing)


    W = Diracoperators.Dirac_operator(U, univ.η, univ.fparam)
    bicg(
        univ.η,
        W',
        univ.φ,
        eps = univ.fparam.eps,
        maxsteps = univ.fparam.MaxCGstep,
        verbose = univ.kind_of_verboselevel,
    )
end

function construct_fermionfield_η!(
    univ::Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp},
) where {Gauge,Lie,Fermi<:StaggeredFermion,GaugeP,FermiP,Gauge_temp} # eta = Wdag^{-1} phi
    U, _... = calc_smearingU(univ.U, univ.fparam.smearing)


    if univ.fparam.Nf == 4 || univ.fparam.Nf == 8
        W = Diracoperators.Dirac_operator(U, univ.η, univ.fparam)
        bicg(
            univ.η,
            W',
            univ.φ,
            eps = univ.fparam.eps,
            maxsteps = univ.fparam.MaxCGstep,
            verbose = univ.kind_of_verboselevel,
        )
    else
        WdagW = Diracoperators.DdagD_operator(U, univ.η, univ.fparam)
        N = get_order(univ.fparam.rhmc_action)

        x = univ.η
        vec_x = univ._temporal_fermi[end-N+1:end]
        for j = 1:N
            Fermionfields.clear!(vec_x[j])
        end

        # eta =  (MdagM)^{-alpha/2} phi -> phi = (MdagM)^{alpha/2} eta
        #(MdagM)^{alpha/2} eta ~ α0 eta + sum_k αk (MdagM + βk)^{-1} eta
        vec_β = get_β_inverse(univ.fparam.rhmc_action)
        vec_α = get_α_inverse(univ.fparam.rhmc_action)
        α0 = get_α0_inverse(univ.fparam.rhmc_action)
        shiftedcg(
            vec_x,
            vec_β,
            x,
            WdagW,
            univ.φ,
            eps = univ.fparam.eps,
            maxsteps = univ.fparam.MaxCGstep,
        )
        Fermionfields.clear!(univ.η)
        Fermionfields.add!(univ.η, α0, univ.φ)
        for j = 1:N
            αk = vec_α[j]
            Fermionfields.add!(univ.η, αk, vec_x[j])
        end
    end
end





end
